Envío del PR inicial: 24.08.2018
Aceptación del PR: 02.09.2018
Escriban una función newton que, a partir de f, fprime y x0 dados, obtenga una de las raices de la ecuación. Comprueben que funciona con $f(x)=x^2 -2$ y $f'(x) = 2x$. Tengan suficiente cuidado para que no haya ningún tipo de inestabilidad de tipo en su función.
Documenta la función de manera adecuada (docstrings).
¿Cómo se comporta, en términos del número de iterados, la convergencia del método de Newton?
Para implementar el método unidimensional, se puede hacer lo siguiente:
"""
newton(f, fprime, x0, número_iteraciones)
`newton` es una implementación unidimensional compleja del método de Newton para encontrar raíces de la función `f`.
# Argumentos
Para poder utilizar la función se requieren los siguientes argumentos:
* `f`, la función compleja (real) de variable compleja (real) de la que se quiere buscar una raíz (un punto \$ a \$ tales que \$ f(a) = 0 \$),
* `f_prime`, la función derivada de `f`, y,
* `x0`, una adivinanza inicial sobre la posición de la raíz.
Opcionalmente, se puede especificar lo siguiente:
* `número_iteraciones`, el número de iteraciones a realizar. (Por defecto está configurado en 1000 iteraciones.)
`newton` requiere que tanto `f` como `fprime` sean funciones, que `x0` sea un número y que `número_iteraciones` sea un número entero. `x0`, en particular, siempre es convertida a un número flotante para mejorar la estabilidad de tipo.
# Ejemplos
```julia-repl
julia> newton(x -> x^2 - 2, x -> 2*x, 2, 5)
1.4142135623730951
julia> newton(x -> x^2 - 2, x -> 2*x, 1.4, 100)
1.414213562373095
julia> newton(x -> x^2 - 2, x -> 2*x, φ, 1)
1.4270509831248424
julia> newton(x -> x^2 - 2, x -> 2*x, φ, 10)
1.414213562373095
julia> newton(z -> z^2 + 1, z -> 2*z, -0.1*im, 5)
0.0 - 1.0032578510960606im
julia> newton(z -> z^2 + 1, z -> 2*z, -0.1*im, 10)
0.0 - 1.0im
```
"""
function newton(f::Function, fprime::Function, x0::Number, número_iteraciones::Int = 1000)
punto_actual = float(x0)
for i in 1:número_iteraciones
punto_actual -= f(punto_actual)/fprime(punto_actual)
end
return(punto_actual)
end
Considerando las funciones $f(x) = x^ 2 - a$ y $g(x) = f'(x) = 2x$ se tiene que si $a \geq 0$, el método de Newton puede encontrar aproximaciones a las raíces cuadradas de $a$ según la adivinanza inicial otorgada.
#Generando el ejemplo con los siguientes parámetros:
f(x) = x^2 - 2
g(x) = 2*x
x0 = 1.5
iteraciones = 10
#La aproximación a la raíz cuadrada positiva de 2 calculada con el método es:
@show raíz_newton = newton(f, g, x0, iteraciones)
#Que el número calculado sea una aproximación a la raíz positiva es resultado de la elección de x0.
#La aproximación a la raíz cuadrada (positiva) de 2 calculada con la función `sqrt` es:
@show raíz = sqrt(2)
#Comparando si son iguales como números flotantes de 64 bits:
@show raíz_newton == raíz
Por lo que el método implementado calculó correctamente la raíz cuadrada positiva de 2 como número flotante de 64 bits con 10 iteraciones, con la adivinanza inicial igual a 1.5.
Para ver la estabilidad de tipo en un ejemplo, basta con usar el macro @code_warntype en la llamada a la función:
#Para el análisis se consideran las funciones y el número de iteraciones definidas para el ejemplo anterior. Lo que se va a modificar va a ser el tipo de adivinanza inicial suministrada.
x0 = 1.5 #Lo mismo que en el ejemplo anterior.
@show typeof(x0) #Número flotante de 64 bits.
@code_warntype newton(f, g, x0, iteraciones)
#El macro no emitió advertencias. Intentando con un entero:
x0 = 2
@show typeof(x0) #Número entero de 64 bits.
@code_warntype newton(f, g, x0, iteraciones)
#El macro tampoco emitió advertencias. Intentando con un irracional:
x0 = MathConstants.φ #El número áureo
@show typeof(x0) #Número irracional
@code_warntype newton(f, g, x0, iteraciones)
#Tampoco emitió advertencias, esto debido al comportamiento de un irracional ante la función float(): se toma la representación como número flotante de 64 bits más cercana a dicho irracional.
float(MathConstants.φ)
# Finalmente, intentando con un complejo de punto flotante:
x0 = 2.0im
@show typeof(x0) #Número complejo flotante de 64 bits
@code_warntype newton(f, g, x0, iteraciones)
Tampoco arrojó advertencias. Estos ejemplos otorgan evidencia sobre la estabilidad de tipo de la función implementada. En todos los ejemplos se devolvió un número flotante de 64 bits.
Para analizar la convergencia del método, consideremos el ejemplo en el que se calcula una aproximación a la raíz cuadrada de dos con una adivinanza inicial constante. ¿Cómo cambia la diferencia entre el valor dado por Julia y el valor calculado con el método de Newton según el número de iteraciones?
Para hacer esto de forma automática y así poder graficar, conviene definir:
Una función que otorgue la diferencia entre el valor de la raíz cuadrada de un número $a$ calculada por el método de Newton implementado y el valor de la misma raíz cuadrada calculada por Julia con una adivinanza inicial fija y un número de iteraciones dado.
Un conjunto de valores de esta diferencia para diferentes cantidades de iteraciones.
"""
diferencia_raíz_cuadrada(a, x0, número_iteraciones)
`diferencia_raíz cuadrada` calcula el valor absoluto de la diferencia entre el valor de la raíz cuadrada positiva del número real positivo \$ a \$ calculada mediante la función `sqrt` implementada en Julia y y la aproximación al valor de la misma raíz cuadrada calculada mediante el método de Newton con el número de iteraciones dado.
# Argumentos
Para poder utilizar la función se requieren los siguientes argumentos:
* `a`, el número real positivo del que se está calculando la raíz,
* `x0`, una adivinaza de la raíz (positiva) de \$ a \$, y,
* `número_iteraciones`, el número de iteraciones realizadas en el método de Newton para calcular la raíz.
#Ejemplos:
```julia-repl
julia> diferencia_raíz_cuadrada(2, 2, 1)
0.08578643762690485
julia> diferencia_raíz_cuadrada(2, 2, 3)
2.1239014147411694e-6
julia> diferencia_raíz_cuadrada(π, 1.5, 1)
0.02474370029108175
julia> diferencia_raíz_cuadrada(π, 1.5, 3)
8.183899780078718e-9
```
"""
function diferencia_raíz_cuadrada(a::Real, x0::Real, número_iteraciones::Int)
raíz_cuadrada_newton = newton(x -> x^2 - a, x -> 2*x, x0, número_iteraciones)
raíz_julia = sqrt(a)
diferencia = raíz_julia - raíz_cuadrada_newton
valor_absoluto = abs(diferencia)
return(valor_absoluto)
end
#El conjunto de números de iteraciones para el cálculo muestra será de la forma:
conjunto_iteraciones = 1:10
#Generando el conjunto de valores para analizar lo que sucede. Se ocupa:
x0 = 100
#Esto debido a que se quiere evitar que la diferencia sea cero para poder usar una gráfica semilogarítimica para visualizar los datos:
valores_raíz_newton = [diferencia_raíz_cuadrada(2, x0, i) for i in conjunto_iteraciones]
#Se requiere cargar la siguiente paquetería:
using Plots
scatter(valores_raíz_newton, title = "Convergencia del método a sqrt(2)", label = "x0 = $x0", yscale = :log10, xlabel = "Número de iteraciones", ylabel = "Diferencia entre valores")